Nonmonotonic critical temperature in superconductor/ferromagnet bilayers 
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The critical temperature T c ol a superconductor/ferromagnet (SF) bilayer can exhibit nonmono- 
tonic dependence on the thickness df of the F layer. SF systems have been studied for a long time; 
according to the experimental situation, the "dirty" limit is often considered which implies that 
the mean free path in the layers is the second smallest spatial scale after the Fermi wavelength. 
However, all calculations reported for the dirty limit were done with some additional assumptions, 
which can be violated in actual experiments. Therefore, we develop a general method (to be exact, 
two independent methods) for investigating T c as a function of the bilayer's parameters in the dirty 
case. Comparing our theory with experiment, we obtain good agreement. In the general case, we 
observe three characteristic types of T c (df) behavior: 1) nonmonotonic decay of T c to a finite value 
exhibiting a minimum at particular df, 2) reentrant behavior, characterized by vanishing of T c in a 
certain interval of df and finite values otherwise, 3) monotonic decay of T c and vanishing at finite df. 
Qualitatively, the nonmonotonic behavior of T c (df) is explained by the interference of quasiparticles 
in the F layer, which can be either constructive or destructive depending on the value of df. 

PACS numbers: 74.50.+r, 74.80.Dm, 75.30.Et 



I. INTRODUCTION 

Superconductivity and ferromagnetism are two com- 
peting orders: while the former "prefers" an antiparallel 
spin orientation of electrons in Cooper pairs, the latter 
forces the spins to align in parallel. Therefore, their co- 
existence in one and the same material is possible only 
in a narrow interval of parameters; hence the interplay 
between superconductivity and ferromagnetism is most 
conveniently studied when the two interactions are spa- 
tially separated. In this case the coexistence of the two 
orders is due to the proximity effect. Recently, much 
attention has been paid to properties of hybrid prox- 
imity systems containing superconductors (S) and fer- 
romagnets (F); new physical phenmneiia were observed 
and predicted in these systems.Era'BQHo One of the most 
striking effects in SF layered structures is highly non- 
monotonic dependence of their critical temperature T c 
on the thickness df of the ferromagnetic layers. Exper- 
iments exploring this nonmonotonic behavior were pec-. 
formedj-previouslM on SF multilayers such as Nb/Gd,LJ 
Nb/Fejfv/V-FeJ] and Pb/Fe,E3 but the results (and, in 
particular, the comparison between the experiments and 
theories) were not conclusive. 

To perform reliable experimental measurements of 
T c (df), it is essential to have df large compared to the 
interatomic distance; this situation can be achieved only 
in the limit of weak ferromagnets. Active experimen- 
tal investigations of SF bilayers and multilayers based on 
Cu-Ni dilute ferromagnetic alloys are carried out by sev- 
eral groups OE3 In SF bilayers, they observed nonmono- 
tonic dependence T c (df). While the reasoitjfor this effect 
in multilayers can be the 0-ir transition^ in a bilayer 
system with a single superconductor this mechanism is 
irrelevant, and the cause of the effect is interference of 



quasiparticle, specific to SF structures. 

In the present paper, motivated by the experiments of 
Refs. [ll|Jl^ we theoretically study the critical tempera- 
ture of SF bilayers. Previous theoretical investigations of 
T c in SF structures were concentrated on systems with 
thin or thick layers (compared to the corresponding co- 
herence lengths) ; with SF boundaries having very low or 
very high transparencies; the exchange energy was often 
assumed to be much larger than the critical tempera- 
ture; in addition, the methads for ., solvin g the problem 
were usually approximate.@a'Bll2(ll3llJllalla The parame- 
ters of the experiments of Refs. |j"j] , |l2| do not correspond 
to any of the above limiting cases. In the present pa- 
per we develop two approaches giving the opportunity 
to investigate not only the limiting cases of parameters 
but also the intermediate region. Using our methods, we 
find different types of nonmonotonic behavior of T c as a 
function of df, such as minimum of T c and even reentrant 
superconductivity. Comparison of our theoretical predic- 
tions with the experimental data shows good agreement. 

A number of methods can be used for calculating T c . 
When the critical temperature of the structure is close to 
the critical temperature T cs of the superconductor with- 
out the ferromagnetic layer, the Ginzburg-Landau (GL) 
theory applies. However, T c of SF bilayers may signif- 
icantly deviate from T cs , therefore we choose a more 
general theory valid at— arbitrary temperature — the 
quasiclassical approach.tlKEj Near T c the quasiclassi- 
cal equations become linear. In the literature the emerg- 
ing problem is often treated with the help of the so-called 
"single-mode" approximation,DoEj£j which is argued to 
be qualitatively reasonable in a wide region of parame- 
ters. However, this method is justified only in a specific 
region of parameters which we find below. Moreover, be- 
low we show examples when this method fails even qual- 
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FIG. 1: SF bilayer. The F and S layers occupy the regions 
— df < x < and < x < d a , respectively. 



itatively. Thus there is need for an exact solution of the 
linearized quasiclassical equations. The limiting case of 
perfect boundaries and large exchange energy was treated 
by Radovic et aln 

Based on the progress achieved for calculation of T c 
in SN systems (where N denotes a nonmagnetic normal 
material), EJ we develop a generalization of the single- 
mode approximation — the multi-mode method. Al- 
though this method seems to be exact, it is subtle to 
justify it rigorously. Therefore we develop yet another 
approach (this time mathematically rigorous), which we 
call "the method of fundaiapntaL soliition" . The models 
considered previousljfi'QBEEiljO'tjO correspond to lim- 
iting cases of our theory. A part of our results was briefly 
reported in Ref. |l[ 

The paper is organized as follows. In Sec. II we formu- 
late the Usadel equations and the corresponding bound- 
ary conditions. Section III is devoted to the exact multi- 
mode method for solving the general equations. An al- 
ternative exact method, the method of fundamental so- 
lution, is presented in Sec. fv[ In Sec. ^ we describe 
results of our methods. In Sec7]Vl|, a qualitative explana- 
tion of our results is presented, applicability of the results 
to multilayered structures is discussed, and the use of a 
complex diffusion constant is c ommented upon. Conclu- 
sions are presented in Sec. VII. Appendixes [a|,[b| contain 
analytical results for limiting cases. Finally, technical 
details of the calculations are given in Appendix |c[ 



II. MODEL 

We assume that the dirty-limit conditions are ful- 
filled, and calculate the critical temperature of the bi- 
layer within the framework of the linearized Usadel equa- 
tions for the S and F layers (the domain < x < d s is 
occupied by the S metal , —df < x < — by the F 
metal, see Fig. II]). Near T c the normal Green function is 



G = sgnw„, and the Usadel equations for the anomalous 
function F take the form 



d 2 F s 
dx 2 



0, < x < d s ; 



Aln- 



T 



- (|w„| + iE cx sgnuj n )Ff = 0, 
-df < x < 0; 



A 

|Wnl 



(1) 

(2) 
(3) 



(the pairing potential A is nonzero only in the S part). 
Here £ s = yjD s /2nT cs , = y / D f /2nT cs are the co- 
herence lengths, while the diffusion constants can be ex- 
pressed via the Fermi velocity and the mean free path: 
D = vl/3; LJ n = 7rT(2n + 1) with n = 0,±1,±2, ... 
are the Matsubara frequencies; E cyi is the exchange en- 
ergy; and T cs is the critical temperature of the S material. 
F s if\ denotes the function F in the S(F) region. We use 
the system of units in which Planck's and Boltzmann's 
constants equal unity, h = fc^ = 1 . 

Equations (0)-@ must be supplemented with the 
boundary conditions at the outer surfaces of the bilayer: 



dF s (d s ) _ dF f (-d f ) 



dx 



dx 



0, 



as well as at the SF boundary: 



dF s (0) = 
dx 

dFf(0) 



70 



dF f (0) 
dx 



C/76 



dx 



F,(p)-Ff(p), 7fc 



Ps€s 

Pftf' 

R b A 



(4) 

(5) 
(6) 



Here p s , Pf are the normal-state resistivities of the S and 
F metals, Rb is the resistance of the SF boundary, and A 
is its area. The Usadel equation in the F layer is readily 
solved: 



Ff = C(u n ) cosh(fc/[x + df]) , 
k f - 



(7) 




iE ex sgnu 7 , 



and the boundary condition at x = can be written in 
closed form with respect to F s : 



dFs(0) 
dx 
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76 + Bf{uj n ) 



Fs(0), 



(8) 



Bf = [kf£f ta,nh(kfdf)]' 



This boundary condition is complex. In order to 
rewrite it in a real form, we do the usual trick and go 
over to the functions 



F ± = F(u n )±F(-w n ). 



(9) 



According to the Usadel equations (|l|)-(||), there is the 
symmetry F(—cu n ) = F*(u n ) which implies that F + is 
real while F~ is a purely imaginary function. 
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The symmetric properties of F + and F~ with respect 
to uj n are trivial, so we shall treat only positive u) n . The 
self-consistency equation is expressed only via the sym- 
metric function Ft'- 



u n >0 v 



(10) 



and the problem of determining T c can be formulated 
in a closed form with respect to F s + as follows. The 
Usadel equation for the antisymmetric function F~ docs 
not contain A, hence it can be solved analytically. After 
that we exclude Ff from boundary condition m) and 
arrive at the effective boundary conditions for Fj 7 : 



^=%)f s + (0), 



dx 
where 

W(u n ) = 7 



dF+(d s ) 
dx 



0. 



A s {^ h + ReBf) + 7 



A s \j b + B f | 2 + 7(76 + Re B f 
A s = k s ^ s tanh(fc s d s ) 



(11) 



(12) 



k s — ~r~ 

4s 



1lT r . 



The self-consistency equation ( |Io| ) and boundary condi- 
tions (pr[)-(12), together with the Usadel equation for 
F+: 



d 2 F+ 
dx 2 



^f: 



2A = 



(13) 



will be used below for finding the critical temperature of 
the bilayer. 

The problem can be solved analytically only in lim- 
iting cases (see Appendix |X]). In the general case, one 
should use a numerical method, and below we propose 
two methods for solving the problem exactly. 



III. MULTI-MODE METHOD 



then the self-consistency Eq. ( )l0| ) takes the form (5 and 
do not depend on u> n ) 



where tj) is the digamma function. 

Boundary condition (fLlT) at x — yields 



51 tan £1 



W(w n ). 



(17) 



(18) 



The critical temperature T c is determined by Eqs. 

Although this method is popular, it is often used with- 
out pointing out the limits of its applicability. We present 
the explicit formulation of the corresponding condition: 
the single-mode method is correct only if the parame- 
ters are such that W can be considered w„-independent 
[because the [left-hand side of Eq. ( |l8| ) must be w„- 
independentlMl 

Appendix g demonstrates examples of the SMA valid- 
ity and corresponding analytical results. 

In one of experimentally relevant cases, E ey _/iTT ca > 1, 
d f ~ £/, the SMA is applicable if ^/E cx /-kT cs > 1/tj, 
(see Appendix |B] for details) . 



B. Inclusion of other modes 

The single-mode approximation implies that one takes 
the (only) real root of Eq. (p~7|) . An exact (multi- 
mode) method for solving problem (|Io|)-(^3|) is obtained 
if we also take imaginaoz-roots into account — there is 
infinite number of these E3 

Thus we seek the solution in the form 

F+(x,w n ) = / (w„)cos 



rn — 1 



fm{Wn)- 



cosh 



cosh f^mff-J 



(19) 



A. Starting point: the single-mode approximation 
and its applicability 

In the single-mode approximation (SMA) one seeks the 
solution of the problem (^0|)-(p^|) in the form 



(14) 
(15) 



Fj{x,U> n ) = /(W„)C0S O 



. x — d. 



A(x) = 8 cos I O 



This anzatz automatically satisfies boundary condition 
(0) at x = d s . 

The Usadel equation (p~2|) yields 



25 

to n + Q 2 irT cs ' 



(16) 



A(x) = Sq cos Qq 



m— 1 



cosh ( 



c — d s 



cosh 



(20) 



(The normalizing denominators in the cosh-terms have 
been introduced in order to increase accuracy of numer- 
ical calculations.) This anzatz automatically satisfies 
boundary condition (11) at x = d s . 

Substituting the anzatz [Eqs. ([l9j)-(|oj)] into the Us- 
adel equation (|l3|), we obtain 

/oK) = (21) 



fm(w n ) — 



- Q 2 TTT r .. 



m = 1,2,. 



4 



then the parameters fl are determined by the self- 
consistency equation ([l0|) (S and £1 do not depend on 

tOn)'. 



T 



i n 2 T r: 



(22) 



From Eqs. ( p2| ) and properties of the digamma function^ 
it follows that the parameters f2 belong to the following 
intervals: 



o < n 2 < 



1 



(23) 



^(2m-l) <n 2 m < ^(2m + l), m=l,2,. 



where 7 E ~ 1.78 is Euler's constant. 

Boundary condition (|ll|) at x — yields the following 
equation for the amplitudes 6: 



So 



W(uj n ) cos (flpds/^s) - ilg sin 

W(uj n ) + VL m tanh (Cl m d s /^ s ) 



m— 1 



0. (24) 



The critical temperature T c is determined by Eqs. ( |22| ) 
and the condition that Eq. (^J) has a nontrivial (w n - 
independent) solution with respect to d. 

Numerically, we take a finite number of modes: m = 
0, 1, . . . , M. To take account of ^-independence of the 
solution, we write down Eq. (|24|) at the Matsubara fre- 
quencies up to the iVth frequency: n — 0, 1, . . . , N. Thus 
we arrive at the matrix equation K nm 5 m = with the 
following matrix K: 



K n = 



u> n /irT cs + fi§ 
W(u> n ) + Q m tanh (fl m ri s 
uj n /nT cs - VL 2 m 
n = 0,1,... ,N, m= 1,2,... ,M. 



(25) 



We take M = N, then the condition that Eq. (|24|) has a 
nontrivial solution takes the form 



det K = 0. 



(26) 



Thus the critical temperature T c is determined as the 
largest solution of Eqs. (E2),(Bh). 



( |l3| ) satisfies the same equations, but with the delta- 
functional "source" O 



„ a _ d 2 G(x,y) 



dx 2 



uj n G(x,y) = -5(x - y), 



(27) 



^^ = ^K) G (0, y) , ^=0. (28) 



(/.?: 



(i:r 



The fundamental solution can be expressed via solutions 
Vx, V2 of Eq. (27) without the delta- function, satisfying 
the boundary conditions at x = and x — d s , respec- 
tively: 



G(x,y;uj n ) = — 



k s /uj n 



smh(k s d s ) + (W/k s S, s ) cosh (k s d s ) 

< Mx \ V2{ y } ; x<y , (29) 

V2{x)vi(y), ys^x 



where 

vi(x) — cosh(k s x) + (W/ks^s) sinh(fc s x), (30a) 
V2 {x) = cosh {k s [x — d s ] ) . (30b) 

Having found G(x, y;uj n ), we can write the solution of 
Eqs. P)-(pl) as 



F+(x;w n ) 



G{x 1 y;u n )A(y)dy. 



(31) 



Substituting this into the self-consistency equation fllOj), 
we obtain 



A(x) In 



T, 



2irT c Y, 



UJ n >0 



A(x) 



G{x,y\w n )A(y)dy 



(32) 



This equation can be expressed in an operator form: 
Aln(T cs /T c ) = LA. Then the condition that Eq. |^) 
has a nontrivial solution with respect to A is expressed 
by the equation 



det L - 1 In — 



0. 



(33) 



The critical temperature T c is determined as the largest 
solution of this equation. 

Numerically, we put problem (|32]),(^3|) on a spatial 
grid, so that the linear operator L becomes a finite ma- 
trix. 



NUMERICAL RESULTS 



IV. METHOD OF FUNDAMENTAL SOLUTION 

By definition, the fundamental solution G(x,y;uj n ) 
(which is also called the Green function) of problem (p"l])- 



In Sees. Ill, IV we developed two methods for cal- 
culating the critical temperature of a SF bilayer. Spec- 
ifying parameters of the bilayer we can find the criti- 
cal temperature numerically. It can be checked that the 
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(nm) 

FIG. 2: Theoretical fit to the experimental data of Ref. [y]. In 
the experiment, Nb was the superconductor (with d a — 11 nm, 
T cs = 7K) and Cuo.43Nio.57 was the weak ferromagnet. From 
our fit we estimate E ex » 130 K and 75 « 0.3. 



multi-mode method and the method of fundamental solu- 
tion yield equivalent results. However, at small temper- 
atures T c <C T cs , the calculation time for the multi-mode 
method increases. Indeed, the size of the matrix K [Eq. 

is determined by the number N of the maximum 
Matsubara frequency ujn, which must be much larger 
than the characteristic energy ttT cs ; hence N ^> T cs /T c . 
Therefore, at low temperatures we use the method of 
fundamental solution. 



A. Comparison with experiment 

Using our methods we fit the experimental data of 
Ref. [n]; the result is presented in Fig. g. Estimating 
the parameters d s = 11 nm, T cs = 7K, p s — 7.5^Slcm, 
£ s = 8.9 nm, p/^j 60/if2cm, £f = 7.6 nm, 7 = 0.15 from 
the experiment and fitting only E cyi and 75, we find 
good agreement between our theoretical predictions and 
the experimental data. 

The fitting procedure was the following: first, we de- 
termine £ cx w 130 K from the position of the minimum of 
T c (df); second, we find 7^ ~ 0.3 from fitting the vertical 
position of the curve. 

The deviation of our curve from the experimental 
points is small; it is most pronounced in the region of 
small df corresponding to the initial decrease of T c . This 
is not unexpected because, when df is of the order of a 
few nanometers, the thickness of the F film may vary sig- 
nificantly along the film (which is not taken into account 
in our theory) , and the thinnest films can even be formed 
by an array of islands rather than by continuous mate- 
rial. At the same time, we emphasize that the minimum 
of T c takes place at df « 5 nm, when with good accuracy 
the F layer has uniform thickness. 




4 d f / A ex 

FIG. 3: Characteristic types of T c (df) behavior. The thick- 
ness of the F la yer is measured in units of the wavelength A cx 
defined in Eq. (M). The curves correspond to different values 
of 7b. The exchange energy is _E ox = 150 K; the other param- 
eters are the same as in Fig. |^. One can distinguish three 
characteristic types of T c (df) behavior: 1) nonmonotonic de- 
cay to a finite T c with a minimum at particular df ("fb = 2; 0.5; 
0.1; 0.07), 2) reentrant behavior (j b = 0.05; 0.02), 3) mono- 
tonic decay to T c = at finite df (7;, = 0). The bold points 
indicate the choice of parameter corresponding to Fig. n. 



B. Various types of T c (df) behavior 

The experimental results discussed above represent 
only one possible type of T c {df) behavior. Now we ad- 
dress the general case; we obtain different kinds of T c (df) 
curves depending on parameters of the bilayer. 

To illustrate, in Fig. |we plot several curves for various 
values of 7?, [we recall that 7^ oc where Rb is the 
resistance of the SF interface in the normal state — see 
Eq. (g)]. The exchange energy is E cx = 150 K; the other 
parameters are the same as in Fig. |[ 

We observe three characteristic types of T c (df) behav- 
ior: 1) at large enough interface resistance, T c decays 
nonmonotonically to a finite value exhibiting a minimum 
at a particular df, 2) at moderate interface resistance, 
T c demonstrates the reentrant behavior: it vanishes in a 
certain interval of df, and is finite otherwise, 3) at low 
enough interface resistance, T c decays monotonically van- 
ishing at finite df. A similar succession of T c (df) curves 
as in Fig. |^ can be obtained by tuning other parameters, 
e.g., the exchange energy E cx or the normal resistances 
of the layers (the parameter 7). 

A common feature seen from Fig. |^ is saturation of 
T c at large df > A cx - This fact has a simple physical 
explanation: the suppression of superconductivity by a 
dirty ferromagnet is only due to the effective F layer with 
thickness on the order of A ex , adjacent to the interface 
(this is the layer explored and "felt" by quasiparticles 
entering from the S side due to the .proximity effect). 

It was shown by Radovic et alii that the order of 
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4 d f l A ex 



FIG. 4: Change of the phase transition's order. This phe- 
nomenon manifests itself as discontinuity of T c (df): the criti- 
cal temperature jumps to zero abruptly without taking inter- 
mediate values. Formally, T c becomes a double-valued func- 
tion, but the smaller solution is physically unstable (dotted 
curve). For illustration we have chosen the curve from Fig. |^ 
corresponding to fb = 0.05. 
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the phase transition may change in short-periodic SF 
superlattices, becoming the first order. We also ob- 
serve this feature in the curves of types 2) and 3) men- 
tioned above. This phenomenon manifests itself as dis- 
continuity of T c (df): the critical temperature jumps to 
zero abruptly without taking intermediate values (see 
Figs. ||,(§). Formally, T c becomes a double- valued func- 
tion, but the smaller solution is physically unstable (dot- 
ted curve in Fig. |]). 

An interesting problem is determination of the tricriti- 
cal point where the order of the phase transition changes. 
The corresponding result for homogeneous bulk super- 
conductors with internal exchange field was obtained a 
long time ago in the framework of the Ginzburg-Landau 
theory.o However, the generalization to the case when 
the GL theory is not valid is a subtle problem which has 
not yet been solved. We note that the equations used in 
Refs. |3|,|l5] were applied beyond their applicability range 
because they are GL results valid only when T c is close 
to T rx . 



FIG. 5: Comparison between single- and multi-mode meth- 
ods. The parameters are the same as in Fig. |j| Generally 
speaking, the results of the single-mode and multi-mode (ex- 
act) methods are quantitatively and even qualitatively dif- 
ferent: b), c), d), and e). However, sometimes the results 
are close: a) and f). Thus the single-mode approximation 
can be used for quick estimates, but reliable results should 
be obtained by one of the exact (multi-mode or fundamental- 
solution) techniques. 



rameters. Using the methods developed in Sees. [in.IV, 
we can check the actual accuracy of the single-mode ap- 
proximation. The results are presented in Fig. |5|. 

We conclude that although at some parameters the re- 
sults of the single-mode and multi-mode (exact) methods 
are close (Figs. || a,f), in the general case they are quan- 
titatively and even qualitatively different [Figs. ||b,c,d,e 
- these cases correspond to the most nontrivial T c (df) 
behavior] . Thus to obtain reliable results one should use 
one of the exact (multi-mode or fundamental-solution) 
techniques. 



Comparison between single- and multi-mode 
methods 



A popular method widely used in the literature for 
calculating the critical temperature of SF bi- and multi- 
layers is the single-mode approxima tion. The condition 
of its validity was formulated in Sec. Ill A. However, this 



approximation is often used for arbitrary system's pa- 



D. Spatial dependence of the order parameter 

The proximity effect in the SF bilayer is characterized 
by the spatial behavior or the order parameter, which 
can be chosen as 



F{x,T = 0) = Tj2F{x,Lo n ), 



(34) 
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FIG. 6: Spatial dependence of the order parameter normal- 
ized by its value at the outer surface of the S layer. Two cases 
are shown differing by the thickness of the F layer df (and by 
the corresponding T c ) at 75 = 0.05. The other parameters are 
the same as in Fig. pi where the chosen cases are indicated 
by the bold points. Although the critical temperatures differ 
by more than the order of magnitude, the normalized order 
parameters are very close to each other, which means that the 
value of T c has almost no effect on the shape of F(x,t — 0). 
The jump at the SF interface is due to its finite resistance. 
With an increase of df the order parameter starts to oscillate, 
changing its sign (this can be seen for the dotted curve, al- 
though negative values of the order parameter have very small 
amplitudes) . 



where r denotes the imaginary time [in the S metal 
F(x,t — 0) oc A(x)]. This function is real due to the 
symmetry relation F(— oj n ) — F*(uj n ). 



We illustrate this dependence in Fig. [| which shows 
two cases differing by the thickness of the F layer df (and 
by the corresponding T c ). Although the critical temper- 
atures differ by more than the order of magnitude, the 
normalized order parameters are very close to each other, 
which means that the value of T c has almost no effect on 
the shape of F(x, r = 0) . Details of the calculation are 
presented in Appendix 



Another feature seen from Fig. g is that the order pa- 
rameter in the F layer changes its sign when the thickness 
of the F layer increases (this feature can be seen for the 
dotted curve, although negative values of the order pa- 
rameter have very small amplitudes). We discuss this 
oscillating behavior in the next section. 
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FIG. 7: Four types of trajectories contributing (in the sense 
of Feynman's path integral) to the anomalous wave function 
of correlated quasiparticles in the ferromagnetic region. The 
solid lines correspond to electrons, the dashed lines — to 
holes; the arrows indicate the direction of the velocity. 



VI. DISCUSSION 



Qualitative explanation of the nonmonotonic 
T c (df) behavior 



The thickness of the F layer at which the minimum of 
T c (df) occurs, can be estimated from qualitative argu- 
ments based on the interference of quasiparticles in the 
ferromagnet. 

Let us consider a point x inside the F layer. Accord-, 
ing to Feynman's interpretation of quantum mechanics ,EZI 
the quasiparticle wave function may be represented as a 
sum of wave amplitudes over all classical trajectories; the 
wave amplitude for a given trajectory is equal to exp(iS'), 
where S is the classical action along this trajectory. We 
are interested in an anomalous wave function of corre- 
lated quasiparticles, which characterizes superconductiv- 
ity; this function is equivalent to the anomalous Green 
function F{x). To obtain this wave function we must 
sum over trajectories that (i) start and end at the point 
x, (ii) change the type of the quasiparticle (i.e., convert 
an electron into a hole, or vice versa). There are four 
kinds of trajectories that should be taken into account 
(see Fig. [?]). Two of them (denoted 1 and 2) start in the 
direction toward the SF interface (as an electron and as 
a hole), experience the Andreev reflection, and return to 
the point x. The other two trajectories (denoted 3 and 
4) start in the direction away from the interface, experi- 
ence normal reflection at the outer surface of the F layer, 
move toward the SF interface, experience the Andreev 
reflection there, and finally return to the point x. The 
main contribution is given by the trajectories normal to 
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the interface. The corresponding actions are 

S% = —Qx — a, (35) 

5*2 = Qx — a, (36) 

S 3 = -Q(2d f + x) - a, (37) 

S 4 = Q(2d f + x) ~ a (38) 

(note that x < 0), where Q is the difference between 
the wave numbers of the electron and the hole, and 
a = arccos(_E/A) is the phase of the Andreev reflec- 
tion. To make our arguments more clear, we assume 
that the ferromagnet is strong, the SF interface is ideal, 
and consider the clean limit first: then Q — k e — — 
y/2m(E + E m + fj,) - y/2m(-E - E cx + JT) « 2E e */v, 
where E is the quasiparticle energy, \x is the Fermi en- 
ergy, and v is the Fermi velocity. Thus the anomalous 
wave function of the quasiparticles is 



ric multilayered structures in the 0-state such as SFS and 
FSF trilayers, SFIFS and FSISF systems (I denotes an ar- 
bitrary potential barrier), and SF superlattices. In such 
systems an SF bilayer can be considered as a unit cell, 
and joining together the solutions of the Usadel equa- 
tions in each bilayer we obtain the solution for the whole 
system (for more details see Sec. VIII of Ref. p8| ). 

Our methods can be generalized to take account of pos- 
sible superconductive and/or magnetic 7r-states (when A 
and/or E ex may change their signs from layer to layer). 
In this case the system cannot be equivalently separated 
into a set of bilayers. Mathematically, this means that 
the solutions of the Usadel eq uations lose their purel y 
cosine form [see Eqs. (@), @, ©, ©, (fog )] 

acquiring a sine part as well. 



F(x) oc exp(iS n ) oc cos(Qdf) cos (Q[df + x]) . (39) 

71 = 1 

The suppression of T c by the ferromagnet is determined 
by the value of the wave function at the SF interface: 
F(0) oc cos 2 (Qdf). The minimum of T c corresponds to 
the minimum value of F(0) which is achieved at df = 
■k/2Q. In the dirty limit the above expression for Q is 
replaced by 



2vr 

Xe-r 



(40) 



(here we have defined the wavelength of the oscillations 
A ex ); hence the minimum of T c (df) takes place at 



,(min) _ 7T 
% f - 9 



El 

E C x 



A, 



ex 

T' 



(41) 



For the bilayer of Rcf. [ll] we obtain df s» 7nm, 
whereas the experimental value is 5 nm (Fig. |^) ; thus 
our qualitative estimate is reasonable. 

The arguments given above seem to yield not only the 
minimum but rather a succession of minima and maxima. 
However, numerically we obtain either a single minimum 
or a minimum followed by a weak maximum (Fig. |^). 
The reason for this is that actually the anomalous wave 
function not only oscillates in the ferromagnetic layer but 
also decays exponentially, which makes the amplitude of 
the subsequent oscillations almost invisible. 

Finally, we note that our arguments concerning oscil- 
lations of F(x) also apply to a half-infinite ferromagnet, 
where we should take into account only the trajectories 1 
and 2 (see Fig. 0). This yields F(x) oc cos(Qa;) (another 
qualitative explanation of this result can be found, for 
example, in Ref. H4) . 



B. Multilayered structures 

The methods developed and the results obtained in 
this paper apply directly to more complicated symmet- 



C. Complex diffusion constant? 

Finally, we comment on Refs. ^|, ^5|Jl^ , |29| , where the 
authors considered (in the vicinity of T c ) diffusion equa- 
tions with a complex diffusion constant D f for the F part 
of the structure. This implies small complex corrections 
to Df over E cx t < 1 in the Usadel equations (r is the 
time of the mean free path). However, we disagree with 
this method for the following reason: although the com- 
plex Df can indeed be formally obtained in the course of 
the standard derivation of the Usadel equationaiJ from 
the Eilenberger onesEEl by expanding over the spherical 
harmonics, one can check that the higher harmonics ne- 
glected in the derivation have the same order of magni- 
tude as the retained complex correction to Df. Hence 
the complexity of Df in the context of the Usadel equa- 
tions is the excess of accuracy. Below we present our 
arguments. 

We give a brief derivation of the Usadel equations 
showing how the complex diffusion constant can be ob- 
tained and why this result cannot be trusted. In the 
"quasi-one-dimensional" geometry (which means that 
the parameters vary only as a function of x) the linearized 
Eilenberger equation in the presence of disorder and the 
exchange field has the form 



>cos6> d 
~~2 dx' 



-F 



Tr +lE -) F = A + ^ (42) 



where, for simplicity, we assume a positive Matsubara 
frequency u> n > 0, and 8 is the angle between the x axis 
and the direction of the Fermi velocity v, while (...) 
denotes angular averaging over the spherical angles. The 
disorder is characterized by the time of the mean free 
path t and the mean free path I (to be used below). In 
the dirty limit the anomalous Green function F is nearly 
isotropic. However, to obtain the Usadel equation for the 
isotropic part of F, we must also take into account the 
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next term from the full Legendre polynomial expansion: 



VII. CONCLUSIONS 



F(x,u n , ( 



k=0 



F k (x,uj n )Pk(cos9) 



F (x,uj n ) + Fi(x,L0 n ) cost 



(43) 



Here we have neglected the harmonics with k ^ 2 assum- 
ing them small; we shall check this assumption later. 

Averaging Eq. ([l||) over the spherical angles first di- 
rectly and second — after being multiplied by cos#, we 
arrive at 



-— Fi 

6 dx 

v d 
2dx~ F ° 

Equation ( f45|) yields 
Fx = 



iE ex ) F = A, 



f 

I 



iE r 



s 1 + 2ljJ n T 

then Eq. @ leads to 
~2~dx 



F x = 0. 



-F , 
dx 



2 F ~(Lu n + iE cx )F Q + A 



0. 



(44) 
(45) 

(46) 

(47) 



D 



vl/3 



1 + 2u)„ T - 



2iE cx r 



Now we must check that the assumption \F\/Fq\ <C 1, 
IF2/J1I <C 1, etc. that we used is indeed satisfied. From 
Eq. (EQ) we obtain 



l/L 



max(l, 2w„r, 2E ox t) 



(48) 



where L is the characteristic space scale on which Fq 
varies. According to the Usadel equation (47), it is given 
by 



(49) 



v/ max(l, 2uo n T, 2E cx t) max(2o-> rl r, 2E cx t) ' 



and the condition of the Usadel equation's validity is 
written as 



/ max(2w n r, 2E ex r) 
max(l, 2w„r, 2E cx t) 



< 1 



(50) 



[similarly, we can also keep the term with k = 2 in series 
(||), which yields \F2jF 1 \ ~ \F X /F \, etc.]. 
Finally, condition (I5Q) takes the form 



2ttT cs t < 1, 



2E CX T < f 



(51) 



(we have taken into account that the characteristic en- 
ergy is oj n ~ irT cs ). 

Now we can analyze our results. If condition (|5l]) is 
satisfied and the Usadel equation is valid, the neglected 
angular harmonics have the relative order of magnitude 
I-F2/-P0I ~ max(27r T cs t, 2E cx t); hence we cannot retain 
the terms of the same order in the diffusion constant [see 
Eq. (p7|)], and we should use the standard expression 
D = vl/3. 



In the present paper we have developed two methods 
for calculating the critical temperature of a SF bilayer 
as a function of its parameters (the thicknesses and ma- 
terial parameters of the layers, the quality of the inter- 
face) . The multi-mode method is a generalization of the 
corresponding approach developed in Ref. 20 for SN sys- 
tems. However, the rigorous justification of this method 
is not clear. Therefore, we propose yet another approach 
— the method of fundamental solution, which is mathe- 
matically rigorous. The results demonstrate that the two 
methods are equivalent; however, at low temperatures 
(compared to T cs ) the accuracy requirements are stricter 
for the multi-mode method, and the method of funda- 
mental solution is preferable. Comparing our method 
with experiment we obtain good agreement. 

In the general case, we observe three characteristic 
types of T c {df) behavior: 1) nonmonotonic decay of T c 
to a finite value exhibiting a minimum at particular df, 

2) reentrant behavior, characterized by vanishing of T c 
in a certain interval of df and finite values otherwise, 

3) monotonic decay of T c and vanishing at finite df. 
Qualitatively, the nonmonotonic behavior olT c {df) is ex- 
plained by interference of quasiparticles in the F layer, 
which can be either constructive or destructive depend- 
ing on the value of df. 

Using the developed methods we have checked the ac- 
curacy of the widely used single-mode approximation. 
We conclude that although at some parameters the re- 
sults of the single-mode and exact methods are close, in 
the general case they are quantitatively and even qual- 
itatively different. Thus, to obtain reliable results one 
should use one of the exact (multi-mode or fundamental- 
solution) techniques. 

The spatial dependence of the order parameter (at the 
transition point) is shown to be almost insensitive to the 
value of T c . 

The methods developed and the results obtained in 
this paper apply directly to more complicated symmet- 
ric multilayered structures in the 0-state such as SFS 
and FSF trilayers, SFIFS and FSISF systems, and SF 
superlattices. Our methods can be generalized to take 
account of possible superconductive and/or magnetic 7r- 
states (when A and/or E ox may change their signs from 
layer to layer). 

We argue that the use of the complex diffusion constant 
in the Usadel equation is the excess of accuracy. 

In several limiting cases, T c is considered analytically. 
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APPENDIX A: ANALYTICAL RESULTS FOR A 
THIN S LAYER 

(i) When d s < £ s and E cx > irT cs , problem @ 
(O) can be solved analytically. The first of the above 
conditions implies that A can be considered constant, 
and F + weakly depends on the spatial coordinate; so 
F + (x,Lu n ) = 2A/uJn+A(ui n ) cosh(k s [x— d s }). The bound- 
ary conditions determine the coefficient A; as a result 



F+(uj n ) = F+(x = 0,uj n ) = — 



A s (u n ) 



A s (co n ) + W(uj n ) 



(Al) 



where k s , A s , and W are defined in Eq. (|12|). Finally, 
the self-consistency equation for T c takes the form 



In- 



F. 



Reip 



1,76 



1 



Tr. 



2 2d slb + B f T c 



(A2) 



where Bf does not depend on oj n due to the condition 
E ex 3> 7rT cs : 



B f = [kf^ftaah(k f df) 



1 



ttTc. 



(A3) 



(i i) If the F layer is also thin, df <C \J D f /2E CX , Eq. 
( |A2| ) is further simplified: 



T r . V \2 To 



-i + T f E c 



E cx 

2-KF r 



"^2 



where r s , Tf are defined similarly to Ref. |28| : 
2d s R b A 2d f R b A 



PsF>s 



Pf D f 



(A4) 



(A5) 



and have the physical meaning of the escape time from 
the corresponding layer. They are related to the quanti- 
ties 7, 76 used in the body of the paper as 



7b 1 d s 
7 7rT cs £ s 



1 d 



Tf = 7b 



(A6) 



(iii) If the S layer is thin, d s <C £ s , and the SF interface 
is opaque, 75 — > 00, the critical temperature of the bilayer 



only slightly deviates from F cs . In this limit Eq. (Al) 
applies with W = "f/"fb ^ 1, & n d we finally obtain: 



T P . = Tc. 



TV 

4t7 



(A7) 



Interestingly, characteristics of the F layer (df, E cx , etc.) 
do not enter the formula.- In particular, this formula is 
valid for an SN bilayertirEj (where N is a nonmagnetic 
normal material, E cx = 0) because Eq. ( |A7j ) was ob- 
tained without any assumptions about the value of the 
exchange energy. 



1. Transparent interface 

When both layers are very thin [d s <C 

with uj r 



y/D s /2w D , 

df <C mm(y/Df/2u) D , y / Df/2E cx ), with uj d the Debye 
energy of the S metal] and the interface is transparent, 
the bilayer is equivalent to a homogeneous superconduct- 
ing layer with internal exchange field. This layer is de- 
scribed by effective parameters: the pairing potential 
A( cff ), the exchange field E^\ and the pairing constant 
\( eS ). In this subsection we develop the ideas of Ref. 53, 
demonstrate a simple derivation of this description, and 
find the limits of its applicability. 

The Usadel equations (|l|),(j|) for the two layers can be 
written as a single equation: 

D f 9(-x) + D s 9(x) d 2 F 
2 ~da? 
- \uo n \F - iE cx sgn(u n )6(-x)F + A8(x) = 0, (A8) 

where 6 is the Heaviside function [6(x > 0) = 1, 0(x < 
0) = 0]. The self-consistency equation (^) can be rewrit- 
ten as 



Ala;) x, U) n ), 



(A9) 



where A is the pairing constant. 

First, we consider the ideal SF interface: 75 = [see 
Eq. (H)], then F(x) is continuous at the interface and 
nearly constant across the whole bilayer, i.e., F s (x\ 



Ff(x) — F. Applying the integral operator to Eq. (|A8|) 



s d s + Vfdf 



dx 



v,d. 



'fdf Jo 



dx (A10) 



(here v is the normal- metal density of states), and can- 
celling gradient terms due to the boundary condition 
we obtain the equations describing a homogeneous layer: 

-\u n \F(u n ) - iEgV sgn(uj n )F(uj n ) + A^ = 0, 

(All) 

A (cff) = x^ cS KtJ2 f M, (A12) 
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with the effective parameters (see also Ref. 



A (eff) = 



-E„ 



A(° ff ) 



(A13) 



-A, ri° ff ) = 



exp 



1 



where 7 E is Euler's constant and Tjf is the critical tem- 
perature of the layer in the absence of ferromagnetism 



(i.e., at E, 



(off) 



0). The critical temperature is deter- 



mined by the equation 



In- 



T, 



(eff) 



fi 4f r 



4> 



(A14) 



Actually, the description in terms of effective parame- 
ters (A13) is applicable at an arbitrary temperature (i.e., 
when the Usadcl equations arc nonlinear) and has a clear 
physical interpretation: the superconducting (A, A) and 
ferromagnetic (E ex ) parameters are renormalized accord- 
ing to the part of time spent by quasiparticles in the 
corresponding layer. This physical picture is based on 
interpretation of r as escape times, which we present in 
the next subsection. 

Now we discuss the applicability of the above de- 
scription for a nonideal interface (76 7^ 0). In this 
case F is nearly constant in each layer, but these con- 
stants arc different: F s (x) ~ F s + C s (x — d s ) 2 , Ff(x) « 
F f +Cf(x+d f ) 2 , where \F S \ > \C s \d 2 s and \F f \ > \C f \dj. 
Using the Usadel equation ( |AS| ) and the boundary con- 
ditions (||),@, we find the difference 6F = F s — Ff. 



5F 



A 



1 

T, 



1 



Tf (\U1„\+ iE ex Sg!lU! n ) 



(A15) 



Finally, the 
\5F/F\ < 1 
yields: 



homogeneous description is val id when 
with F determined by Eq. (All)], which 



max(B cx , lu d ) max(r s , Tf) < 1 



(A16) 



(here uj n ~ u) D has been taken as the largest characteris- 
tic energy scale in the quasi-homogeneous bilayer). 



2. Interpretation of r as escape times 

The quantities r s , Tf introduced in Eq. (^5|) may be 
interpreted as escape times from the corresponding lay- 
ers. The arguments go as follows. 

If the layers are thin, then the diffusion inside the lay- 
ers is "fast" and the escape time from a layer is deter- 
mined by the interface resistance. The time of penetra- 
tion through a layer or the interface is determined by the 
corresponding resistance: R s (f) or Rf,, hence the diffusion 
is "fast" if < R b . 

Let us use the detailed balance approach, and consider 
an interval of energy dE. In the S layer, the charge in 



this interval is Q s = ev s dEAd s . Let us define the escape 
time from the S layer t s , so that the current from S to F 
is equal to Q s /t s . On the other hand, this current can 
be written as dE/eRb, hence 



and we immediately obtain 



dE 
ei? 6 ' 



t» = 



d s R b A 
p s D s 



(A17) 



(A18) 



Similarly, we obtain the expression for the escape time 
from the F layer tf. As a res ult, the relations between 
the quantities r defined in Eq. (A5) and the escape times 
t are simply 



2t, 



r f = 2tf. 



(A19) 



Microscopic expressions for the escape times may be 
obtained using the Sharvin formula for the interface re- 
sistance. Assuming, for definiteness, that the Fermi ve- 
locity is smaller in the S metal, v s < Vf, we obtain 



R b = 



e 2 v s v s X 



and consequently 



t s = tt — r b , 



tf 



V f d f 



(A20) 



(A21) 



where r b is the inverse transparency of one channel. The 
asymmetry in these expressions stems from our assump- 
tion v s < V f . In t he opposite case the indices s and / in 
Eqs. flA20D,flA2"l|) should be interchanged. 



APPENDIX B: APPLICABILITY OF THE 
SINGLE-MODE APPROXIMATION 



As pointed out in Sec. Ill A , the single- mode approxi- 
mation (SMA) is applicable only if the parameters are 
such that W [see Eq. (|l2|)] can be considered ui n - 
independent. An example is the case when 75 ^> |-B/|, 
hence W = 7/76. 

The condition 7^ 3> \Bf\ can be written in a simpler 
form; to this end we should estimate \Bf\. We intro- 
duce the real and imaginary parts of kf. kf = k'j + ik'L 
and note that k'f > k'j. Then using the properties of 
the trigonometric functions and the estimate tanhx ~ 
min(l,a;) we obtain 



\B f 



[k'fCft^h(k'fdf)] \ 



(Bl) 



and finally cast the condition 7^ 3> \Bt\ into the form 



lb 



<C mm ■ 



/ max 



T r Er. 



T ' 7rT 







T r E r 



(B2) 
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where the ratio T c /T cs originates from uj n /wT cs with 
io n ~ 7rT c as the characteristic energy scale in the bi- 
layer. 



If condition (B2) is satisfied, then the SMA is valid 
and T c is determined by the equations 



T r ^ \ 2 2 T r 



o, ai) ( n^-) = 1 

W 7b 



(B3) 
(B4) 



These equations can be further simplified in two lim- 
iting cases which we consider below. 

(1) -Lf « 1: 



in this case Eq. ((bJ) yields ft 2 = and Eq. @ 

takes the form 



hi- 



T 

J- cs 



4> 



2 2 7h d s T c 



(B5) 



which reproduces the 7;, ^> \Bf\ limit of Eq. (A2) 



7b 5s 

in this case Eq. (g) yields f2|j = §, and Eq. @ 



takes the form 



2 8 



2 rt i 2 



T 

-L CS 

~Tr. 



(B6) 



Equations (B3)-(B6) can be used for calculating the 
critical temperature T c and the critical thickness of the S 

(cr) 

layer d s below which the superconductivity in the SF 
bilayer vanishes (i.e., T c = 0). 



1. Results for the critical temperature 



In the limit when T c is close to T cs , Eqs. flB5D ,flB6D 
yield 



T c = T„ [ 1 
and 



4 7b d t 



■t 7 „ • (ds ^ 
it — <C mm 



7b 



6 <^s 



(B7) 



1- 



4 i 



if ds ^ (i 
it — 3> max 1, — 

6 V 7 



(B8) 



Using relations ( |A6|) one can check that result (B7) is 
equivalent to Eq. (A.7). 



2. Results for the critical thickness 



(cr) 

The critical thickness of the S layer d s is defined as 
the thickness below which there is no superconductivity 



in the SF bilayer: T c (di cr) ) = 0. When T c -> 0, Eq. 



yields f2 = 1/ A /27 F (where 7 B as 1.78 is Euler's 
constant), and Eq. (B4) takes the form 



tan 



(cr) 

Explicit results for d s can be obtained in limiting cases: 



1 dj ciy 



1_ 
lb ' 



(B9) 



6 7b 7b Cs 



and 



7 (cr) 



if^»l. 

7b 6 



(BIO) 



(Bll) 



APPENDIX C: SPATIAL DEPENDENCE OF THE 
ORDER PARAMETER 

According to the self-consistency equation, in the S 
layer the order parameter F (x, t — 0) is proportional to 
A(x): 



F s {x,t = 0) = 



AW 
ttA 



(CI) 



where A is the pairing constant which can be expressed 
via the Debye energy: 



ln 



( 2 1e^d 
V 7TT 0S 



(C2) 



The pairing potential A(x) can be found as the eigenvec- 
tor of the matrix L — l\n(T cs /T c ) [see Eq. (0)], corre- 
sponding to the zero eigenvalue. 

After that we can express F(x, t = 0) in the F layer 
via A(x) in the superconductor. The Green function 
Ff(x,u) n ) in the F layer is given by Eq. (ffl), with C(oj n ) 
found from the boundary conditions: 



B 



J 



F s (0,oj n ) 



(C3) 



v 7b + Bf J cosh(kfdf) 
The Green function at the S side of the SF interface is 
F+(0,u n )+F-(0,LJ n ) 



F s (0,uj n 



(C4) 



The symmetric part F^~ is given by Eq. (pl|) . The anti- 
symmetric part is 

F- = C~(u; n )cosh(k s [x ~ d s }) , (C5) 

with C~(uj n ) found from the boundary conditions: 



47 Im B f 



A s |7 b + %| 2 + 7 ( 7b + Re%) 



F+(0,O _ 
cosh(k s d s ) 
(C6) 
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Finally, the order parameter in the F layer is the Fourier 
transform [see Eq. (fsi])] of 



Ff{x,u. 



1 



i"f Im Bf 



A a \<y b + B / \* + >y(rn ) + ILeB f ) 
B f \ cosh(k f [x + d f ]) 



Ib + Bj 



Mkfdf) 



G(0,y;cj n )A(y)dy. 

(C7) 
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